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Abstract 

This article consists of a brief discussion of the energy density over time or frequency 
that is obtained with the wavelet transform. Also an efficient algorithm is suggested 
to calculate the continuous transform with the Morlet wavelet. 

The energy values of the Wavelet transform are compared with the power spec- 
trum of the Fourier transform. Useful definitions for power spectra are given. 

The focus of the work is on simple measures to evaluate the transform with the 
Morlet wavelet in an efficient way. The use of the transform and the defined values 
is shown in some examples. 
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1 Introduction 



The wavelet transform is a method for time-frequency analysis. The applica- 
tion for acoustic signals can be found in several publications. These publica- 
tions deal for example with the analysis of dispersive waves [Tf2] . source or 
damage localization [3f4"] . investigation of system parameters [5116] or active 
control [7] . A comparison of the short time Fourier transform and the wavelet 
transform is done by Kim et.al. [8]. It is found that the continuous wavelet 
transform CWT of acoustic signals is a promising method to obtain the time 
- frequency energy distribution of a signal. 

Another sort are wavelets for spatial transforms, that are an effective tool for 
damage localisation [9lfT0] . 

The aim of this publication is to present an algorithm for the continuous 
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wavelet transform CWT with the Morlet wavelet. Caprioli et al. state in [TT] 
that " The high computational complexity (of the CWT) is not a serious worry 
for today's computer power; still it might be an obstacle for an "on line" ap- 
plication." Alternatives to the CWT are characterised by low computational 
complexity, but the best results are obtained by the CWT. The presented al- 
gorithm implements simple measures to reduce the computational complexity 
of the CWT. An implementation in Java is written and is open source and 
freely available onlinj^l. 

Nevertheless it should be mentioned that the Wigner approximation is also 
very useful in the given context. A recent publication is for example [12j . Other 
methods, that are adapted to the analysis of dispersive waves are discussed in 

pirn]. 

This publication starts with a discussion of the energy values that are obtained 
by the wavelet transform. It follows a description of an efficient algorithm to 
evaluate the transform. 



2 Morlet Wavelet Transform 



The theoretical background of the wavelet transform can be found in textbooks 
[T5][T6] . Only the used definitions are stated. The wavelet transform with the 
wavelet ip of a signal y{t) is defined as 



WJ(0, b) = -j== f yQtW (—) dt, (1) 

where a is called dilatation and b translation parameter. The Morlet wavelet 
■0 which is sometimes called Gabor wavelet is given by 

xj;(t) = e-^e juJot (2) 



and c$ = yir/(3. The values (3 = and uj are defined in a particular ap- 
plication so that the admissibility condition is valid [16]. The function (J2J) is 
called mother wavelet. 

The Morlet wavelet is common for the time frequency analysis of acoustic sig- 
nals. For a comparison of different so called mother wavelets see Schukin et. 
al. IE 
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2.1 Comparing Fourier and wavelet transform 



When using the Fourier transform one usually develops the spectrogram. The 
wavelet transform has scaling factors. The analog to the spectrogram is the 
scalogram defined as 

m(a,b)\ 2 . (3) 
The scalogram is a measure of the energy distribution over time shift b and 
scaling factor a of the signal. It holds that the energy E y of a signal y is 

oo oo oo 

E y = J \y(t)fdt= J J \W^aM^F- (4) 

— oo — oo — oo 

If instead of the scaling factor a the frequency value / = 1/a is used, the value 
/ is only the real frequency if uj = 2n. 
It follows with da = jjf df = — J? df that 

oo oo 

E y = J J \Wl{fM 2 dfdb. (5) 

— oo — oo 

It is possible to divide this total energy into an energy density over time and 
over frequency. This is achieved by one integration over frequency or time. 
The energy density over time is defined by 

oo 

Et(b)= J \W»(f,b)\ 2 df. (6) 

— oo 

The energy density over frequency, or the energy density spectrum is defined 
by 



E f (f)= \WlUM Z db. (7) 



2.1.1 Power signals 

The value \W^(f, b)\ 2 is sometimes called Morlet power spectrum, but here 
the term energy density is used. A Comparision of the values obtained with 
the Fourier transform and the wavelet transform is, for example, given in |18j . 
A power signal is characterized by 

1 °° 

P = _ J \ x {t)\ 2 dt < oo. (8) 

— oo 

If the wavelet transform is applied to a power signal one recognises that for 
higher frequencies the energy density is lower, as shown in the examples in 
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section HI To achieve a better match, Shyu [19] proposes a modified equal 
amplitude wavelet power spectrum that is given by 

MPS = ^|^(/,6)| 2 . (9) 

The above definition corresponds more closely to power spectra obtained with 
the discrete Fourier transform. This can be explained since P = and c^a 
is the effective length of the wavelet so T e jf = c^/f is the value to scale the 
energy density. Since one is usually familiar with power spectra, the above 
definition (Q is useful but a problem is that 

oo oo 

J J MPS dfdb. (10) 

— oo — oo 

This fact can be compared with the windowed Fourier transformation. If the 
window like for example, the flat top is defined so that the peak value is 
const antf^l. it follows that the effective noise band width is altered by the 
windowing function. The Morlet wavelet transform can be interpreted as a 
windowed Fourier transformation with a frequency dependent window size. 



2.1.1.1 Alternative definition of the power spectrum In the follow- 
ing a new definition of the power spectrum is given. It is defined analogous to 
the discrete Fourier transform, so that the power can be calculated by sum- 
mation of the power values. The discrete Fourier transform is a filter with the 
bandwidth Af, so each value of the power spectrum is the integral over Af. 
The power spectrum that is defined with the energy values of the discrete 
continuous Morlet wavelet transform is 

Pi = ^Ef(fi)- (11) 

Given the shape of the Morlet wavelet, the only sensible scaling of the fre- 
quency values is logarithmic. With this prerequisite the values follow the 
same tendency as equation flU]). Using equation (TTTT) the total power is given 
by P - ^ P>. 

From equation fTTTl) and At = T/N t , where N t is the number of time values, 
it follows that 

Pd(i,j) = ^ Afi At \W^(fi, bj)\ 2 = ^\WlU u b^ (12) 

is a reasonable definition for the power spectrum over time. An example for a 
power signal is given in section 14.11 

2 Which is strictly true only if f n = T/(2mr), but the window should minimize the 
deviation. 



4 



2.1.2 Energy signals 

However, the discrete Fourier Transform is usually the best method to work 
with power signals. The wavelet transform in contrast is suited to work with 
energy signals, for which 

E = J \x(t)\ 2 dt < oo (13) 

— oo 

is true. It is useful to use the value Ef for comparing Fourier and Wavelet 
spectra in this case. This is done by calculating 

W) = P^ (14) 

from the power spectrum obtained with the Fourier transform. An example of 
an energy signal is given in 14.21 



3 CWT Algorithm 



Usually the continuous wavelet transform is implemented in a direct algorithm. 
The starting point is a program written by Hayashi [20] , where the integral in 
equation ([1]) is implemented by a simple summation over the elements. Typical 
values for the scaling parameter a are 20 to 40, where n , with a = \/2. 

The time shift parameter is b m = mAt, where < m < N . 
The algorithm to calculate the wavelet transform of a discrete signal yi of the 
length N is for each a n : 

(1) Calculate ipi, of the length 2N, b m = mAt with ujq = V2n/ At. 

(2) Shift t/> about m : Ef Vi ' i>N- m +i 

The frequency is given by / = Tp 3 -. 



3.1 Improved algorithm 



The continuous wavelet transform is very inefficient compared with the dis- 
crete wavelet transform. A brief description of a new algorithm follows which 
introduces ideas that are crucial for the high efficiency of the discrete wavelet 
transform DWT. Besides the following improvements the program calculates 
the corresponding frequency values and energy or power distributions that are 
defined in section 12.11 
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3.1.1 Effective Compact support 



The outer parts of the wavelet \t\ ^> have no significant contribution to the 
result because of the exponentially decaying envelope term e _/3t I 2 which acts 
like a window, known as the windowed Fourier transform. Due to the limited 
floating point precision of a computer, the outer parts will not contribute to 
the integral in equation ([1]). 

An a-priori estimation of the effective wavelet length is made here. For t = 
the envelope has its maximum which is unity. The length is chosen so, that for 
all \t\ > \t max \ the envelope is smaller than the precision. For a 64 bit floating 
point number this means that 



e -&?nax/2 = 2.22 ■ 1(T 16 

1^1=8.49/^/3. (15) 

Strictly the outer parts can still contribute to the integral (CD) especially if y 
is a transient function, but in practise even a weaker definition does not show 
a significant effect. 

The concept of neglecting the outer parts is called effective compact support. 
Note that to develop the transform the wavelet is scaled with t/a. A general 
straight forward implementation is not possible, due to the discrete time val- 
ues. This is a reason why the discrete wavelet transform uses scaling values 
that are defined by a n = 2 n , where n is an integer value. 
The actual implementation works in the following way. The length of the 
wavelet that represents the lowest frequency f m i n is chosen so that T/ At = 2 N . 
The mixed approach is that the wavelet length is halved for each doubling of 
the frequency, so for each value /„ = 2 n f min , and left otherwise unchanged. 
For all frequency values /, for which /„ = 2 n f min < fi < 2 n+1 f min = f n+1 
holds, it follows that the support of the wavelet is even longer than t max . The 
reason for this effect is that also t max is scaled and so for a constant length 
and a higher frequency the support of the mother wavelet is longer. 



3.1.2 Frequency dependent time shift 

Due to the argument — of the wavelet ip in ([T]) the transform can be inter- 
preted as a convolution of a signal y and a wavelet ip 

oo 

y**P= J y{t)i>{t - t) dt. (16) 

— oo 

With the scaled argument t/a and r = b/a equation ([1]) follows. The algorithm 
that is outlined in section [3] does actually not set the translation parameter b, 
but b/a. It is the highest precision of the continuous translation parameter b 
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with discrete time values. But it is inefficient, since the whole wavelet and its 
time resolution are scaled with a. The resolution depends on the wavelength. 
If one increases the time shift for lower frequencies this does not result in a 
less accurate resolution. Like for the effective compact support the translation 
parameter is adjusted for each doubling of the frequency. 
The implementation uses the frequency dependent length of the a period T\ = 
1//. The time shift parameter b is defined relative to the period length T\ and 
an arbitrary dyadic factor 2 Ub in a way that T\j2 Ub > b > T\/2 nb+1 is valid. 
The details can be found in the Java-code. 

In some cases the values given above are not fully applicable. For example for 
frequencies / > f s /2 nb , with the sampling frequency f s , the smallest possible 
shift is At which can be bigger than T\/2 nb , but it is the highest possible 
precision. The default value is 2™ 6 = 16, which results in a high resolution 
where a difference to a higher resolution is not recognizable. With this measure 
a reduced computational effort is achieved. 

Together with the compact support it is possible to use an amount of 200 
scaling factors instead of 20 and wait a few seconds for the transform to be 
calculated. 



3.1.3 Frequency dependent truncated boundaries 

One may interpret y(t) as an infinite signal of which a part of the length T is 
known. The convolution of the wavelet over the signal results in a considerable 
violation of the admissibility condition since all the parts of the wavelet that lie 
outside the known signal must vanish. In practise this results in invalid data at 
the boundaries. When working with the compact support this is avoided if the 
algorithm guaranties that the whole wavelets always fits in the signal. Since 
the compact support of the wavelet will be shorter for high frequencies this 
results in a frequency dependent truncation. This measure is also implemented 
in [21]. 



3.2 Fast Morlet wavelet transform 



A fast implementation is described in [2T] and often called fast Morlet wavelet 
transform. The convolution in equation (jTJ) is implemented by building the 
DFT of the wavelet and the analysed function, multiplying them and build- 
ing the inverse Fourier transform. Since the DFT is calculated with the fast 
Fourier transform (FFT) algorithm, the computational effort is greatly re- 
duced compared with [2D] . 
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3.3 Computational effort 

This section will compare the order of the amount of operations that are 
needed for each implementation. Since the amount of operations depends on 
the actual implementation such an investigation is beyond the scope of this 
investigation. 

AGU-implementation For the implementation [20] the signal of length N 
has to be multiplied N times with the wavelet. For M scaling factors the 
order of the amount of operations is 0(M x N 2 ). 

Fast Morlet Wavelet transform The FFT has a computational effort of 
0(log(N)N), if iV is a power of two. It follows that the order of operation 
for the fast Morlet wavelet transform is 0(M x log(N)N). 

New algorithm According to section l3.1. li the number of points of the wavelet 
N w is independent of the signal length N. For a given frequency value fi 
the number of points is inversely dependent on the frequency N w ~ 1/ fi. 
According to section l3.1.2l the translations parameter is also inversely depen- 
dent on the frequency foj ~ 1/ fi- The wavelet has to be multiplied N/bi times 
with the signal. Since N w < N the amount of multiplications is N w N/bi. It 
follows that the order of the amount of multiplications for a given frequency 
value is O(N). The resulting overall effort is 0(M x N). 



3.3.1 Benchmark 

The consumed time for the transform depends on many other factors. It is only 
an indication for the amount of operation that actually have to be performed. 
Important factors that effect the computation time is the systems memory 
and the implementation itself. 

The benchmark compares the standard Morlet wavelet transform (st), the fast 
Morlet wavelet transform (fast) and the described new Algorithm (new) in a 
typical application. 

The standard Morlet wavelet transform is implemented with Matlab's conv() 
function, that implements the convolution with a Matlab's f ilter () function. 
Matlab's profiler shows that 99.3% of the computational time is used by the 
function conv() which is optimized. 

The fast Morlet wavelet transform is implemented by a self-written convolu- 
tion function that uses Matlabs highly optimized fft() function. The fft()- 
function itself chooses the presumable best algorithm for the application. The 
FFT algorithm needs a number of points that are an integer power of two 2 n . 
This can be problematic if one wants to adjust the length of a window so that 
it fits the discrete best case frequencies f n = T/(2wr). This defiancy does not 
apply to the continuous wavelet transform. 

The described algorithm is implemented with Java and called from Matlab. It 
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0.17 
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0.16 


4.83 
0.15 


9.12 
0.14 
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0.14 



Table 1 



Computation time for three different algorithms: standard (st), fast and new Morlet 
wavelet transform and for different sizes. 



is not optimised and it is known that Java programs of this sort usually suffer 
a 50-fold performance degradation [22J compared to native Fortran or similar 
applications. The intention of the implementation is to follow the guidelines 
for software engineering: Build the functionality first and do the optimization 
afterwards. The advantages of Java are that it is platform independent, has 
an easy to use interface to Matlab and no special or commercial libraries are 
needed. Table 13.3.11 shows the computational time for the same M = 200 
frequency values and a maximum frequency f max = / s /8. It is assumed that 
measurements are often taken with a portable computer, so a Samsung X20 
notebook is chosen, equipped with Pentium M CPU at 1.73 GHz and 1GB 
main memory. The parameters for the new transform were chosen such that 
the results are identical, a time shift parameter of 7\/4 and t max = 5.6s. The 
results show a good repeatability and the same tendency on other computers. 
The order of operations that is derived in section [3731 is validated with the nor- 
malised computational time T/N 2 for the standard, T / (Nlog(N)) for the fast 
and T/N for the new algorithm. If the computational time would solely de- 
pend on the amount of operations the normalised computational time should 
be constant. 

When using the standard and fast algorithms the size of the memory limited 
the test to a maximum number of 2 15 points. With the new algorithm it was 
possible to calculate the transform also for 2 18 and 2 19 in a time of 37.25s and 
69.16s. These values also indicate that the estimation of a linear increasing 
number of operations is valid, since Tx ^° — is again 0.14. For a typical sam- 
pling frequency of f s = 48kHz this means that it is possible to analyse around 
10s. 
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4 Examples 



The capabilities of the wavelet transform and the described algorithm are 
demonstrated in the examples below. 



4-1 Power signal 



A typical power signal is the sine function. Consider the example with a chang- 
ing frequency 

f Asin(27rfit) for t < 1/2 

y{t) = \ { ! 1 (17) 

[ A sin(27r/ 2 t) for t > 1/2, 

with At = 2~ 13 , /1 = 80Hz, f 2 = QAOHz and A = 4. The power density 
spectrum equation (TTTT) is obtained with the fast Morlet wavelet transform 
and the presented algorithm for 100 frequency values. The results are identical 
to those of the other algorithms and plotted in figure [H 
The peaks are very broad which is expected and due to the e~PI 2t windowing 
term. The peak at j\ = 80Hz is lower than that at 6A0Hz, which is due 
to the truncation described above and that is plotted in figure [2j The peak 
value does not match the amplitude. This is due to the fact that the total 
power is approximately correct. Since the peak is very broad only the peak 
value or the total can power can be correct. The total power is P = J2i Pi — 
8.1 which corresponds to the correct value of P = 1/2 A 2 = 8. Using the 
equal amplitude wavelet power spectrum that is given in equation §§§ and 
performing an integration over time results in a good approximation of the 
peak value. 

The energy density spectrum is plotted in the bottom part of figure [lj showing 
an unusual representation of the example equation (j!7p . 
A contour plot of the power density spectrum over time is given in figure [21 
Note that the high frequency component f 2 is more accurately localized in the 
time domain. 



4-2 Energy signal 



The tested signal is 



. . sin(a/t) 
y(t) = for t min <t<t r 
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Fig. 2. Contour plot of the power density spectrum over time, equation (fT2j) . of the 
example equation (|17p built with the Morlet wavelet 
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thereby, the parameters are set to a = 10, 1/ At = f s = 2 12 , t min = y a/ (27if max ), 
fmax fs/& and t max .25s. 

The energy is plotted in figure El The solid curve is calculated with the dis- 
crete Fourier transform (DFT) and equation ffT4j) . The Fourier transform of 
function ffl8l) can be obtained by the convolution of the Bessel-function J 
which is the Fourier transform of sm ("/ f ) an d the sinc-function which is the 
Fourier transform of the rectangular window of the size T = t max — t min . 



, , 7i , , — , sin(Tc<j) 
F{v(t)} = 7T^o(2V^) * — (19) 



The small fluctuations in figure [3] correspond to the sinc-function and the high 
and long to the Bessel-function. 

The dashed curve is calculated with the Morlet wavelet transform. This curve 
follows more clearly the theoretical dependence that is plotted as the 

dash dotted curve next to the others. 

A contour plot of the energy density spectrum over time is given in figure HI 
It shows good agreement with the actual frequency of this function which is 
given by uj{t) = a ft 2 . 
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Fig. 4. Contour plot of the energy density spectrum over time of equation flTHJ) 
5 Concluding remarks 

The study shows that it is possible to obtain an efficient transform with an 
algorithm that is restricted to the values that are numerically significant. This 
efficiency results in less computational effort and less memory consumption. 
Nevertheless, for applications that are time critical it will be worth to work 
on a optimized version in a native language and to apply optimisation. Since 
the programme is tested and open source, it will facilitate the programming 
of such an optimised version. 

The examples and definitions of energy and power values should help with the 
interpretation of the values obtained with the wavelet transform. 
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